cornelius.senf@geo.hu-berlin.de

The 4 R Sessions:

  1. Basic concept of programming, data types, vectors and factors

  2. Matrices and data.frames, data import/export, graphics and visualization

  3. Tidy data, data wrangling

  4. Today: Spatial data in R

Raster data in R

Raster data in R

library(raster)

Raster data in R

ras <- raster(ncol = 100, nrow = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)

ras
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

Raster data in R

m <- rnorm(n = ncell(ras), mean = 0, sd = 1)

values(ras) <- m

plot(ras)

Raster data in R

Raster data in R

ras2 <- ras + 10

ras3 <- ras * ras2

ras4 <- ras^2

ras5 <- (ras - cellStats(ras, stat = mean)) / 
                                cellStats(ras, stat = sd)

Raster data in R

Raster data in R

ras.na <- ras

ras.na[ras.na > 0.5] <- NA

Raster data in R

Raster data in R

k <- matrix(c(  0, 0.5, 1,
              0.5, 1.5, 2,
              1.5, 6.0, 3), 
            ncol = 3, byrow = TRUE)
k
##      [,1] [,2] [,3]
## [1,]  0.0  0.5    1
## [2,]  0.5  1.5    2
## [3,]  1.5  6.0    3
ras.reclass <- reclassify(x = ras, rcl = k)

Raster data in R

Raster data in R

ras[10, 20]
##           
## -2.136716
ras[10, 20] <- 10

Raster data in R

Raster data in R

ras.stack <- stack(ras, ras2)

ras.stack
## class       : RasterStack 
## dimensions  : 100, 100, 10000, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## names       :   layer.1,   layer.2 
## min values  : -3.497505,  6.502495 
## max values  :   10.0000,   14.1605

Raster data in R

layer1 <- subset(x = ras.stack, subset = 1)

layer1
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer.1 
## values      : -3.497505, 10  (min, max)

Raster data in R

data <- as.data.frame(ras.stack)

head(data)
##      layer.1   layer.2
## 1 -1.2699512  8.730049
## 2  0.5653630 10.565363
## 3 -0.3784131  9.621587
## 4 -0.6302751  9.369725
## 5 -0.1924782  9.807522
## 6  0.2749651 10.274965

Raster data in R

elevation <- raster("Data/elevation.envi")

elevation
## class       : RasterLayer 
## dimensions  : 709, 1306, 925954  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 16.90833, 27.79167, 44.33333, 50.24167  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/corneliussenf/Dropbox/Teching/Quantitative Methods/Week-13/Data/elevation.envi 
## names       : elevation

Raster data in R

projection(elevation)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

See for projection definitions: http://spatialreference.org/

Raster data in R

plot(elevation)

Raster data in R

Raster data in R

climate <- stack("Data/climate.envi")

nlayers(climate)
## [1] 4
names(climate)
## [1] "X.Tmax."  "X.Tmin."  "X.Tmean." "X.Prec."
names(climate) <- c("Tmax", "Tmin", "Tmean", "Prec")

Export raster data

writeRaster(elevation, 
            filename= "Data/elev.envi",
            overwrite = TRUE)

Vector data

Vector data

library(sp)

Vector data: simple random sample

coords_srs <- data.frame(x = runif(100, 17, 27), 
                         y = runif(100, 45, 50))

pts_srs <- SpatialPoints(coords = coords_srs)
pts_srs
## class       : SpatialPoints 
## features    : 100 
## extent      : 17.02223, 26.78149, 45.0296, 49.96563  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

Vector data: Assign Coordinate system

proj <- projection(elevation)

CRS(proj)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
pts_srs <- SpatialPoints(coords = coords_srs,
                    proj4string = CRS(proj))

Vector data

plot(elevation)
plot(pts_srs, add = TRUE)

Advanced plotting

library(mapview)

mapview(climate)

Advanced plotting

Vector data

df.climate <- extract(climate, pts_srs, df = TRUE)

class(df.climate)
## [1] "data.frame"
head(df.climate)
##   ID     Tmax      Tmin    Tmean     Prec
## 1  1 140.4167 32.583332 86.25000 53.66667
## 2  2 105.2500 23.500000 64.16666 64.50000
## 3  3 117.9167 30.000000 73.66666 68.58334
## 4  4 145.7500 40.000000 92.75000 52.66667
## 5  5  91.5000  9.166667 50.16667 71.58334
## 6  6 116.3333 22.916666 69.41666 64.08334

Vector data

spdf.climate <- extract(climate, pts_srs, sp = TRUE)

spdf.climate
## class       : SpatialPointsDataFrame 
## features    : 100 
## extent      : 17.02223, 26.78149, 45.0296, 49.96563  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 4
## names       :             Tmax,              Tmin,            Tmean,             Prec 
## min values  : 61.9166679382324, -16.4166660308838, 22.5833339691162, 43.4166679382324 
## max values  :           166.75,  67.9166641235352,              117, 84.5833358764648

Vector data

You can acess SpatialPointsDataFrame the same way as data.frame!

head(spdf.climate)
##       Tmax      Tmin    Tmean     Prec
## 1 140.4167 32.583332 86.25000 53.66667
## 2 105.2500 23.500000 64.16666 64.50000
## 3 117.9167 30.000000 73.66666 68.58334
## 4 145.7500 40.000000 92.75000 52.66667
## 5  91.5000  9.166667 50.16667 71.58334
## 6 116.3333 22.916666 69.41666 64.08334
spdf.climate$Tmax[1]
## [1] 140.4167

Reading vector data

library(rgdal)

Reading vector data

pts.sample <- readOGR('Data', 'pts_sample')
## OGR data source with driver: ESRI Shapefile 
## Source: "Data", layer: "pts_sample"
## with 142 features
## It has 1 fields
## Integer64 fields read as strings:  id
pts.sample
## class       : SpatialPointsDataFrame 
## features    : 142 
## extent      : 17.80209, 27.04407, 44.64419, 49.84834  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : id 
## min values  :  0 
## max values  : 99

Writing vector data

writeOGR(spdf.climate, "Data", "pts_climate", driver = "ESRI Shapefile", overwrite_layer = TRUE)